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1  Introduction 


Probit  analysis  is  the  primary  statistical  tool  used  to  analyze  dose-response  data  from  laser  retinal 
damage  threshold  experiments.  Beginning  about  1970,  the  U.S.  Army  Medical  Research  Detachment 
(USARMD)  performed  this  analysis  using  small  programs  developed  in-house  and  based  on  the 
procedure  described  by  Frisch  [1],  Since  1996,  USAMRD  has  been  using  computer  programs  that 
have  evolved  from  the  original  EZ-probit  program  developed  by  the  U.S.  Air  Force  Armstrong 
Laboratory  [2],  The  current  version  of  this  program  is  Probit  version  2. 1.2. 3.  Although  it  has 
undergone  minor  corrections  and  revisions,  Probit  has  existed  essentially  in  its  current  form  since 
1998. 

The  author  has  developed  a  new  program,  ProbitFit,  as  a  replacement  for  Probit.  The  primary 
enhancement  offered  by  ProbitFit  is  an  on-screen  graph  of  the  dose-response  data,  as  well  as  the  re¬ 
sulting  fits.  In  addition,  a  rudimentary  data  editor  is  included,  allowing  the  researcher  the  capability 
to  rapidly  explore  several  “what  if”  scenarios. 

ProbitFit  was  written  using  the  Java  programming  language,  and  therefore  should  run  on  any 
operating  system  for  which  an  appropriate  Java  Runtime  Environment  is  available.  Although  based 
on  entirely  new  code,  ProbitFit  uses  the  same  iterative  fitting  procedure  developed  by  Finney  [3] 
that  is  used  by  Probit.  Therefore,  as  a  primary  design  goal,  ProbitFit  was  required  to  reproduce 
the  results  of  Probit  when  used  to  process  the  same  data  set. 

This  report  first  provides  a  brief  overview  of  probit  analysis.  The  iterative  routine  used  by  ProbitFit 
to  fit  the  dose-response  data  is  then  outlined,  followed  by  a  description  of  the  fiducial  limit  calcula¬ 
tions.  There  is  no  attempt  at  mathematical  rigor;  for  details,  the  reader  is  referred  to  Finney  [3], 


2  Overview  of  Probit  Analysis 


Probit  Analysis  is  a  procedure  used  to  analyze  data  in  which  the  outcome  of  an  experiment  is 
quantal:  there  is  a  response  (yes),  or  there  is  not  a  response  (no).  The  probability  of  producing 
a  response  depends  on  the  input  parameters  to  the  experiment.  In  an  experiment  to  determine 
the  laser-induced  retinal  damage  threshold  at  a  specified  wavelength,  pulse  duration,  and  irradiance 
diameter,  the  input  parameter  is  the  energy  of  a  pulse  from  the  laser.  A  yes  response  corresponds  to 
the  observation  of  a  lesion  in  the  retina;  no  response  means  no  alteration  of  the  retina  is  observed. 


2.1  Dose- Response  Curve 

In  probit  analysis,  dose-response  data  is  assumed 
for  a  given  dose  q,  the  probability  of  a  response  is 

p(x) = -jL=  r 

V2tvo2  J — oo  l  j 

where  x  =  log10(<j).  The  goal  is  to  find  the  values  of  p  and  a  which  best  describe  a  set  of  experimental 
data. 

Figure  1  shows  the  characteristic  sigmoid  shape  of  this  distribution.  The  curve  in  this  figure  illus¬ 
trates  a  normal  distribution  having  a  mean  p  =  1.0  and  standard  deviation  a  =  0.2. 


to  follow  a  log-normal  distribution.  That  is, 


exp 
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Figure  1:  Example  dose-response  curve  with  a  log-normal  distribution  having  a 
mean  fi  —  1.0  and  a  standard  deviation  a  =  0.2.  These  values  correspond  to 
ED50  =  10M  =  10.0  and  ED84/ED50  =  ED50/ED16  =  10cr  =  1.58.  The  gray  lines 
indicate  the  16%  level  (/x  —  <r),  the  50%  level  (/x),  and  the  84%  level  (/x  +  a). 


The  notation  ED50  is  used  to  refer  to  the  the  median  effective  dose,  that  is,  dose  at  which  the 
probability  of  producing  a  result  is  50%  (e.g.  the  laser  pulse  energy  expected  to  produce  a  retinal 
lesion  in  half  of  the  exposures).  It  is  the  ED5o  value  obtained  from  a  laser  exposure  experiment  that 
is  quoted  as  the  damage  threshold.  From  Equation  1,  we  can  see  that  P(x  =  /x)  —  0.5.  Therefore, 
the  ED50  is  related  to  the  mean  of  the  probability  distribution  as  ED50  =  10M.  The  ED50  extracted 
from  the  curve  of  Figure  1  is  ED50  =  101  =  10.0. 

The  steepness  of  the  curve  in  Figure  1  is  related  to  the  standard  deviation  a  of  the  distribution.  For 
exposure  data  exhibiting  a  very  sharp  threshold,  a  will  be  very  small,  and  the  dose-response  curve 
will  approach  a  step  function.  If  the  exposure  data  is  of  poorer  quality,  a  will  be  larger,  and  the 
dose  response  curve  will  have  a  finite  slope. 

For  the  log-dose  value  x84  =  /x  +  <7,  the  probability  of  producing  a  lesion  is  P{x%i)  =  0.84.  The  ED84 
dose  is  given  by  ED84  =  10^+<T .  The  ratio  EDg4/ED5o  =  1017  is  directly  related  to  the  standard 
deviation  of  the  dose-response  distribution.  In  reports  produced  by  USAMRD,  the  ratio  ED84/ED50 
is  referred  to  as  the  “slope  of  the  probit  curve,”  and  is  used  as  a  measure  of  how  tightly  a  damage 
threshold  has  been  determined  from  laser  exposure  data.  Typically,  EDg4/ED50  is  found  to  be  in 
the  range  1.0  to  2.0. 


Note:  At  the  log-dose  value  x\e  —  fi  —  <7,  the  probability  of  producing  a  lesion  is  P(x ie)  =  0.16. 
In  this  case  ED16  =  10“-<7,  and  therefore  the  ratio  ED50/ED16  =  1017  also  gives  a  direct  measure 
of  the  standard  deviation  of  the  dose-response  distribution. 


2.2  Probit  Transformation 

Making  the  change  of  variable 


(*~M) 

a 


(2) 


2 


in  Equation  1  leads  to  the  following  expression  for  the  probability  distribution: 

(a  —  f») 

By  using  the  probit  transformation,  defined  as 


(3) 


a 


(4) 


the  probability  distribution  (Equation  3)  becomes 


w  =  £ 


du 


(5) 


The  probit  value  Y  =  5  corresponds  to  the  mean  log-dose  p.  Increasing  the  log-dose  value  x  by  lcr 
increases  the  probit  value  by  1. 


Note:  The  Normal  Equivalent  Deviate  is  defined  as  b  =  (x  —  p)/cr,  which  takes  on  negative 
values  for  x  <  p.  This  made  data  analysis  more  difficult,  especially  before  the  wide  availability 
of  desktop  computers.  Defining  the  probit  with  an  offset  value  of  5  generally  limits  the  probit 
values  to  positive  numbers.  A  probit  of  Y  =  0  corresponds  to  x  —  p  —  5er.  Prom  Equation  1, 
P(p— 5<t)  =  2.9  x  10~7,  a  probability  value  unlikely  to  be  encountered  in  any  realistic  experiment. 


Figure  2  shows  the  dose-response  curve  of  Figure  1  replotted  on  a  linear  scale  of  probit  values.  The 
(now  non-linear)  percentage  axis  is  drawn  on  the  right  side  of  the  plot.  The  probit  transformation 
(Equation  4)  has  linearized  the  data;  the  probit  value  is  a  linear  function  of  the  log-dose  value: 


Y  =  a  +  bx 

The  slope  and  intercept  are  (from  Equation  4) 


(6) 


a  =  5  —  p/cr  (7) 

b  —  l/a  (8) 

Linearization  makes  it  possible  to  extract  information  from  the  data  using  a  graphical  analysis  [1] . 
After  drawing  a  best-fit  line  by  eye  on  a  plot  of  probit  vs.  log-dose  (Figure  2),  p  and  cr  are  deter¬ 
mined  using  Equations  7  and  8.  However,  the  probit  transformation  also  serves  as  the  basis  for  the 
automated  fitting  procedure  used  by  ProbitFit  (Section  3). 


Note:  The  Probit  program  reports  the  value  b  from  Equation  6  as  the  “slope  of  the  probit 
curve.”  Prom  Section  2.1,  a  =  log10(EDs4/ED5o).  Recall  that  the  ratio  ED84/ED50  is  the 
USAMRD  definition  of  the  “slope”  (Section  2.1).  The  two  definitions  of  the  “slope”  are  related 
through  Equation  8: 


l°gio  (ED  84  /ed50) 
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0.99 


dose 


Figure  2:  The  dose-response  curve  of  Figure  1  replotted  on  a  linear  scale  of  pro¬ 
bit  values  (Equation  4).  The  right-hand  axis  shows  the  corresponding  percentage 
values.  The  gray  lines  indicate  the  16%,  50%,  and  84%  percentage  levels. 


2.3  Binomial  Distribution  and  Exposures  at  a  Single  Dose 

The  log-normal  distribution  of  Equation  1  gives  the  probability  P{x)  of  a  lesion  being  produced  by 
a  single  exposure.  If  n  exposures  are  made  at  the  same  dose,  the  probability  of  producing  r  lesions 
is  determined  from  the  binomial  distribution: 

pr(r|”)  =  ri(^b)ip(ir(1~'p<1,)""'  (10) 

The  mean  and  variance  in  the  expected  number  of  lesions  is 


f  -  nP(x) 

Var(r)  =  nP(x)(l  —  P(  a:)) 


(ID 

(12) 


2.4  Damage  Threshold  Experiments  and  Probit  Analysis 

An  experimental  determination  of  the  damage  threshold  for  laser-induced  retinal  injury  generally 
consists  of  making  exposures  at  a  number  of  pulse  energies  in  the  retinas  of  one  or  more  test  subjects 
(usually  non- human  primates,  or  NHPs).  The  selected  pulse  energies  are  expected  to  bracket  the 
damage  threshold  value.  Multiple  exposures  at  each  pulse  energy  may  be  made.  The  retinal  exposure 
sites  are  then  examined  ophthalmoscopically  (usually  1  hour  or  24  hours  after  exposure).  The 
presence  of  a  minimal  visible  lesion  (MVL)  is  generally  used  to  determine  whether  there  was  or  was 
not  a  response  to  the  particular  laser  pulse. 

Assume  n,  exposures  are  made  at  the  pulse  energy  (dose)  Qi,  and  r*  lesions  are  observed,  for  each 
of  the  i  =  1  ,...,&  pulse  energies  used  in  the  experiment.  From  Equation  11,  the  ratio  pi  =  Ti/rii 
provides  an  estimate  of  the  dose-response  curve  (Section  2.1)  at  the  log-dose  value  a q  =  log^fe)- 


4 


The  experimentally  determined  dose-response  curve,  that  is,  the  plot  of  the  ratios  Pi  versus  the 
log-dose  values  aq,  is  expected  to  exhibit  the  characteristic  sigmoid  shape  of  Figure  1.  The  goal  of 
the  probit  analysis  is  to  determine  the  values  of  the  parameters  p  and  a  such  that  Equation  1  best 
describes  the  data.  The  value  for  ED50  and  its  fiducial  limits  (Section  4)  may  then  be  determined. 


3  ProbitFit  Fitting  Procedure 


ProbitFit  uses  an  iterative  fitting  procedure  to  find  the  slope  and  intercept  of  the  probit  line.  The 
mathematical  derivation  of  the  procedure,  based  on  the  Method  of  Maximum  Likelihood,  is  described 
in  detail  by  Finney  [3] .  It  is  difficult  to  work  directly  with  the  normal  distribution  curve  of  Equation  1 
to  derive  the  maximum  likelihood  equations.  This  is  because,  for  example,  partial  derivatives  of  P(x ) 
with  respect  to  p  and  <7  are  required.  Therefore,  the  derivation  of  the  fitting  procedure  is  based  on 
using  the  probit  transformation  (Equation  4)  to  linearize  the  data. 

The  total  probability  of  observing  rq  lesions  in  rq  trials  at  log-dose  a q  for  each  of  the  i  =  1, . . . ,  k 
doses  (pulse  energies)  is 


k 

P(ri,...,rfc)  =  JJPr(rq|rq) 


(13) 


It  is  assumed  that  the  “true”  dose-response  distribution  maximizes  P.  The  total  probability  P  is 
proportional  to  exp(L),  where  the  likelihood  function  is 

k  k 

L(a,  b )  =  £r<  log(-PTO)  +  £(«i  -  rO  log(l  -  P(Yi))  (14) 

i=l  i= 1 

The  dependence  of  L  on  the  parameters  a  and  b  enters  through  Equations  5  and  6  for  P{Y). 
Maximizing  the  likelihood  function  L  will  maximize  the  total  probability  P.  Doing  this  leads  to  an 
iterative  procedure  which,  for  each  iteration,  is  identical  to  performing  a  weighted  linear  regression 
for  the  parameters  a  and  b. 


3.1  Input 

The  input  needed  for  the  fitting  procedure  is  the  experimentally-determined  dose-response  data: 


Xi 

rii 

n 

Pi  =  n/ni 


log10  of  the  ith  dose  level 

#  of  trials  at  dose  level  aq 

#  lesions  observed  at  dose  level  Xi 
Experimental  estimate  for  P(aq) 


Note,  however,  that  the  input  to  ProbitFit  consists  of  the  raw  dose-lesion  data,  which  lists,  for  each 
trial,  the  dose  (beam  or  pulse  energy)  at  which  the  trial  was  run,  and  whether  a  lesion  was  produced. 
The  program  will  process  the  raw  data  to  calculate  the  aq,  rq,  rq  and  pi  values. 
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3.2  Iterative  Fitting  Procedure 

The  iterative  fitting  procedure  used  by  ProbitFit  can  be  outlined  as  follows: 


1.  Guess  an  initial  value  for  the  slope  6  and  intercept  a.  ProbitFit  starts  with  6  =  0  and  a  =  5. 


2.  Calculate  Y,  =  a  +  bxi  for  each  log-dose  value  xt. 

3.  Use  Equation  5  to  calculate  Pi  =  P(Yi  —  5)  for  each  value  aq. 


4.  Calculate 


Zi  = 


exp 


-(Xi  -  5) 
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5.  Calculate  the  working  probits 


Vi  =  Yi  +  ( Pi  ~  Pi)/Zi 


(15) 


(16) 


6. 


A  weighted  linear  regression  [4]  of  the  working  probits  yl  versus  the  log-doses  a q  is  performed 
to  find  new,  improved  values  for  the  slope  6  and  intercept  a.  The  weight  factor  for  each  dose 
level  is  n;w ,,  where 


Pi(l-Pi) 


(17) 


7.  Go  back  to  step  2,  repeat  until  convergence. 


The  mean  p  —  log10(ED5o)  (Equation  1)  is  calculated  at  the  end  of  each  iteration.  This  is  done  by 
setting  Y  —  5  in  Equation  6.  Convergence  is  considered  to  occur  when  the  change  in  p  from  one 
iteration  to  the  next  is  smaller  than  a  user-controlled  tolerance  level. 


Note:  ProbitFit  actually  produces  a  fit  of  the  form  Y  =  y  +  b(x  —  x),  where  x  and  y  are 
weighted  averages  of  the  log-dose  and  working  probit  values.  This  is  because  6,  x  and  y  are 
used  to  calculate  fiducial  limits  on  ED50.  (The  weighting  factors  in  calculating  the  averages  are 
equal  to  TiiWi,  where  Wi  is  determined  from  Equation  17.) 


3.3  Goodness  of  Fit,  x2 

Although  the  fitting  procedure  produces  a  linear  fit  of  the  form  of  Equation  6,  recall  that  the  goal 
is  to  fit  the  experimental  data  to  the  log-normal  distribution  curve,  Equation  1.  Therefore,  the 
measure  of  the  goodness  of  the  fit  is 


2  _  V-'  ( Pi  ~  P(xi )) 

X  Var(p0 

P(xi)  is  to  be  calculated  using  the  final  fit  result.  The  weight  of  each  point  (dose  level)  is  determined 
from  the  binomial  distribution.  Since  Pi  =  r,/?ii ,  then  from  Equation  12 


=  -EM 

Tlj 


(19) 
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Therefore 


2 


(pi  -  P(xj))2 
~[Ul  P{Xi){\  -  P{Xi)) 


Using  equations  16  and  17,  this  can  be  reduced  to  a  form  useful  for  calculations: 


k 

X2  =  ^2  n'Wi  [Vi  ~Yi) 
i= 1 


(20) 


(21) 


where  yt  is  the  working  probit,  and  Yl  is  calculated  from  the  final  fit.  A  statistical  test  of  the  quality 
of  the  fit  is  then  provided  by  the  x2  distribution  for  (k  —  2)  degrees  of  freedom,  where  k  is  the 
number  of  different  pulse  energies  used  in  the  experiment. 


4  Fiducial  Limits 


The  fiducial  limits  are  an  estimate  of  the  range  of  values  for  the  ED50  which  is  supported  by  the 
statistical  analysis  of  the  data  with  some  specified  level  of  confidence. 


4.1  Uncertainties  in  the  Fit 

The  probit  fitting  procedure  provides  an  equation  to  calculate  the  expected  probit  value  for  a  given 
log-dose  value  x 


Y  =  y  +  b(x  —  x )  (22) 

The  fitting  procedure  also  provides  estimates  Var(y)  and  Var(6)  for  the  variance  in  the  fit  parameters 
y  and  b.  Prom  this,  the  variance  in  the  expected  probit  value  may  be  determined: 


Var(Y)  =  Var(y)  +  (x  —  x)2Var  (£>)  (23) 

Note  that  Var(Y)  increases  as  the  x  moves  away  from  the  mean  value  x  of  the  log-dose  values  used 
in  the  exposure  experiment. 


4.2  Fiducial  Limits  for  the  Probit  Value 

In  the  usual  statistical  interpretation,  experimentally  determined  probit  values  are  expected  to  be 
normally  distributed  about  a  mean  value  Y  (Equation  22)  with  a  standard  deviation  a  =  yVar(Y). 
68%  of  measured  probit  values  are  expected  to  lie  within  one  <7  of  Y.  One  says  that  the  fiducial 
limits  of  Y  are  given  by  Y  ±  ^/Var(Y)  at  a  confidence  level  of  68%. 

More  generally,  the  upper  and  lower  fiducial  limits  to  Y  may  be  calculated  as 
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Figure  3:  Example  of  a  probit  curve,  including  the  95%  confidence  level  fiducial 
curves.  The  horizontal  gray  line  indicates  the  50%  probability  of  response  level.  As 
indicated,  the  position  at  which  this  level  intersects  the  probit  curve  and  the  two 
fiducial  curves  indicates  the  ED5o  and  the  Upper  and  Lower  Fiducial  Limits  for  the 
ED50. 


Yufl  =  Y  +  ty/Vax(y)  (24) 

Ylfl  =Y-  ty/VriX)  (25) 

where  t  is  the  normal  deviate  for  the  desired  confidence  level.  That  is,  find  t  such  that 

_  .  ,  I  f 1  ,  1  9\  ,  „  „  confidence  level  ,  , 

Pnorm{t )  =  ^  J  exP(~2U  )  du  =  0-5  +  - - -  (26) 

The  value  t  —  1.96  corresponds  to  a  95%  confidence  level. 


The  far  right  side  of  Equation  26  can  be  understood  as  follows:  Pnorm  (1.96)  =  0.975,  which 
means  that  the  probability  of  finding  a  probit  value  Y  >  Yufl  is  2.5%.  By  the  symmetry  of 
the  normal  probability  density  function,  the  probability  of  finding  a  probit  value  Y  <  Ylfl  is 
also  2.5%.  Therefore,  the  probit  value  is  expected  to  be  in  the  range  Ylfl  <  Y  <  Yufl  with  a 
probability  of  95%. 


Figure  3  shows  the  results  of  a  probit  analysis  on  a  set  of  experimental  laser  exposure  data.  The  95% 
confidence  level  fiducial  curves  have  been  included.  Note  that  they  have  the  characteristic  bowed 
shape  expected  from  Equations  24,  25  and  23. 
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4.3  Fiducial  Limits  on  the  Dose 


Of  greater  interest  is  the  dose  and  the  fiducial  limits  on  that  dose  that  produces  a  response  with  a 
specified  probability  level.  Of  particular  interest  here  is  the  laser  beam  power  or  pulse  energy  that 
will  produce  a  retinal  lesion  in  50%  of  exposures,  i.e.,  the  ED50. 

The  method  for  locating  the  fiducial  limits  for  the  ED50  is  illustrated  in  Figure  3.  The  horizontal 
line  in  the  figure  indicates  the  50%  probability  of  response  level.  This  level  intersects  the  straight 
probit  curve  at  the  ED50  dose  level  (as  indicated  by  the  dropped  vertical  line  in  the  figure).  The 
dose  at  which  the  50%  level  intersects  the  fiducial  curve  to  the  left  of  the  probit  line  is  the  Lower 
Fiducial  Limit  for  the  ED50,  while  the  dose  at  which  this  level  intersects  the  fiducial  curve  to  the 
right  of  the  probit  line  gives  the  Upper  Fiducial  Limit  for  the  ED50. 

Mathematically,  Equations  24  and  25  are  inverted  to  solve  for  the  values  xjjfl  and  xlfl  such  that 
Yufl  and  Ylfl  are  equal  to  the  probit  value  corresponding  to  the  desired  response  level.  For  the 
ED50,  this  value  is  Y  =  5.  For  the  response  probability  P,  the  fiducial  limits  on  xp  =  log10(EDp) 
are 


g  . 

XUFL(LFL)  =  XP  +  ^  —  (Xp  -X)± 


*>(1-0) 


(l  —  p)Var(y)  +  [xp  —  x)  Var (b) 


9  =  ^2Var(b) 


(27) 

(28) 


The  +(-)  corresponds  to  the  upper  (lower)  limit.  The  upper  and  lower  fiducial  limits  for  the  EDp 
are  then  l0a:[,Fi  and  1(F''FI'. 

Because  of  the  shape  of  the  fiducial  curves  (and  also  because  the  independent  variable  is  the  log10  of 
the  dose),  the  fiducial  limits  will  be  asymmetric  about  EDp.  As  Figure  3  indicates,  the  asymmetry 
increases  the  further  the  probability  level  is  from  50%. 


If  the  heterogeneity  factor  h  =  x2/(&  —  2)  is  large  (i.e.  is  significant  as  measured  by  the  x2 
distribution  for  k  —  2  degrees  of  freedom,  where  k  is  the  number  of  exposure  levels  considered), 
then  the  variances  in  Equation  27  must  be  multiplied  by  h,  and  the  factor  t  is  obtained  from  a 
Student’s  t-distribution.  (See  Finney  [3]  for  details.) 


5  Comparison  of  the  ProbitFit  and  Probit  Programs 


In  their  report,  Cain,  Noojin  and  Manning  [2]  compared  several  procedures  for  performing  a  probit 
analysis  with  the  Finney  method  (as  implemented  in  Probit).  Although  they  noted  some  differences 
in  the  results  obtained  by  different  methods,  these  differences  were  generally  inconsequential  com¬ 
pared  to  experimental  uncertainties.  In  this  section,  we  compare  results  from  ProbitFit  and  Probit 
for  three  sample  data  sets.  Numerical  values  are  deliberately  quoted  to  more  decimal  places  than 
would  normally  be  reported. 


5.1  70  msec  Pulse  Data 

This  data  set,  taken  from  Table  A-II  of  Frisch  [1],  is  for  70  msec  pulsed  exposure  data.  The  data 
has  been  binned  by  the  pulse  energy.  Table  1  lists  results  of  fitting  this  data  set  using  ProbitFit  and 
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Table  1:  Comparison  of  fits  to  70  msec  pulse  data 


Parameter 

ProbitFit 

Probit  2. 1.2. 3 

ED50 

12.400 

12.400 

Upper  FL 

13.751 

13.751 

Lower  FL 

10.848 

10.848 

ED  84 

18.568 

18.568 

UpperFL 

22.280 

22.280 

Lower  FL 

16.536 

16.536 

Slope  of  fit,  b 

5.67 

5.67 

Iterations 

7 

56 

Figure  4:  Dose-response  and  probability  curve  graphs  of  ProbitFit  fits  to  70  msec  pulsed 
data  [1]. 


Probit.  Figure  4  shows  a  graph  of  the  data  set  and  the  resulting  fit  which  has  been  generated  by 
the  ProbitFit  program. 

The  two  programs  have  produced  identical  results,  although  the  ProbitFit  program  required  far 
fewer  iterations.  The  programs  were  run  with  their  default  settings.  In  particular,  the  default 
tolerance  used  to  test  for  convergence  in  the  fitting  procedure  is  far  tighter  in  Probit  (1  part  in 
10~18  for  Probit  versus  1  part  in  10-9  for  ProbitFit).  The  results  obtained  here  indicate  that  the 
default  tolerance  used  by  Probit  is  far  tighter  than  necessary  to  obtain  an  good  fit  to  the  data. 


5.2  3.5  nsec  Pulse  Data  at  540  nm 

This  data  set  was  obtained  for  exposures  at  the  wavelength  540  nsec  for  3.5  nsec  pulses  [5].  In 
this  case,  the  data  was  not  pre-binned  by  pulse  energy  before  analysis.  A  comparison  of  the  results 
obtained  using  ProbitFit  and  Probit  is  presented  in  Table  2.  Once  again,  the  two  program  produce 
the  same  results.  Graphs  of  the  ProbitFit  results  are  shown  in  Figure  5. 
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Table  2:  Comparison  of  fits  to  3.5  nsec  pulse  data  at  540  nm 


Parameter 

ProbitFit 

Probit  2. 1.2. 3 

ED50 

6.218 

6.218 

Upper  FL 

7.208 

7.208 

Lower  FL 

5.503 

5.503 

ED84 

7.990 

7.990 

Upper  FL 

11.083 

11.083 

Slope  of  fit,  b 

9.13 

9.13 

Iterations 

14 

23 

Figure  5:  Dose-response  and  probability  curve  graphs  of  ProbitFit  fits  to  3.5  nsec  pulsed 
data  at  A  =  540  nm  [5]. 
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Table  3:  Comparison  of  fits  1.33  fim  corneal  threshold  dat a 


Parameter 

ProbitFit 

Probit  2. 1.2. 3 

ED50 

233.188 

233.188 

Upper  FL 

78.987 

78.982 

Lower  FL 

327.742 

327.744 

ED§4 

506.640 

506.640 

Upper  FL 

1453.767 

1432.849 

Slope  of  fit,  b 

2.95 

2.95 

Iterations 

10 

56 

Figure  6:  Dose-response  and  probability  curve  graphs  of  ProbitFit  fits  corneal  threshold 
data  at  A  =  1.33  gm  [6], 


5.3  1.33  // m  corneal  threshold  data 

This  data  set  was  recorded  for  a  measurement  of  the  corneal  damage  threshold  for  exposures  at  a 
wavelength  of  1.33  pm  [6].  A  comparison  of  the  ProbitFit  and  Probit  results  is  presented  in  Table  3. 
Figure  6  shows  a  graph  of  the  ProbitFit  results. 

Although  the  ED50  and  EDg4  values  are  identical,  there  is  a  discrepancy  at  the  10~5  level  between 
the  fiducial  limits  calculated  by  the  two  programs.  The  origin  of  this  discrepancy  is  not  clear.  It 
is  unlikely  to  be  related  to  the  convergence  test  for  the  fit  procedure,  since  the  EDp  values  are 
identical.  The  difference  is  considered  insignificant  for  practical  purposes. 


6  Conclusion 


The  ProbitFit  program  to  perform  a  probit  analysis  of  dose-response  data  has  been  developed  by 
the  author  for  USAMRD.  ProbitFit  produces  essentially  identical  results  as  Probit  version  2. 1.2. 3 
when  used  to  analyze  the  same  data  set,  while  adding  the  significant  capability  to  produce  graphs  of 
the  data  and  resulting  fits.  ProbitFit  can  therefore  be  recommended  as  a  new  tool  for  the  analysis 
of  laser  damage  threshold  experiments. 
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A  Selected  Source  Code  Files 

ProbitFit  was  written  using  the  Java  programming  language.  It  was  developed  using  the  Java 
1.5  JDK  from  Sun  Microsystems  [7].  Much  of  the  source  code  deals  with  the  user  interface,  and 
is  not  reproduced  here.  This  appendix  contains  the  source  code  for  the  classes  used  to  store  data, 
implement  the  iterative  fitting  procedure,  and  subsequently  calculate  the  effective  dose  (EDp)  values 
and  associated  fiducial  limits. 


A.l  TrialDataltem.java 

package  ng.probitfit .data; 

/* 

*  Class  to  hold  the  result  of  a  single  trial,  i.e.,  the  dose  used, 

*  and  whether  or  not  there  a  response  was  observed  for  that  dose. 

*/ 

public  class  TrialDataltem  implements  Comparable<TrialDataItem>  { 

public  double  dose; 
public  boolean  response; 

/* 

*  Constructors 
*/ 

public  TrialDataltem  ()  { 
dose  =0.0; 
response  =  false; 

> 

public  TrialDataltem  (double  d,  boolean  r)  { 
dose  =  d; 
response  =  r; 

> 

/* 

*  String  representation  for  printing  to  file 
*/ 

public  String  toString  0  { 

return  String. format  ("7,-10. 3f  y,d",  dose,  response  ?  1  :  0)  ; 

> 

/* 

*  Implementation  of  Comparable  interface 
*/ 

public  int  compareTo  (TrialDataltem  tdi)  { 

//  order  determined  by  ’dose’  field 
return  Double . compare  (dose,  tdi. dose); 

> 

> 


14 


A. 2  TrialData.java 


package  ng.probitfit .data; 
import  java. util.*; 

/* 

*  A  collection  of  TrialDataltem  objects 
*/ 

public  class  TrialData  extends  ArrayList<TrialDataItem>  { 

//  default  data  description 

private  static  String  DEF_DESCRIPTION  = 

//  description  of  the  data 

private  String  description  =  DEF_DESCRIPTION; 

/* 

*  Clear  the  data 
*/ 

public  void  clear  ()  -[ 
super. clear  (); 

description  =  DEF_DESCRIPTION; 

> 

/*  Set  or  get  the  description  of  the  data  */ 
public  void  setDescription  (String  desc)  { 
description  =  desc; 

> 

public  String  getDescription  ()  -[ 
return  description; 

> 

public  void  sort  ()  { 

Collections . sort  (this); 

> 

> 


A. 3  DoseResponseDataltem.java 

package  ng.probitfit .data; 

/* 

*  Class  to  hold  a  single  probit  data  point  of  the  dose-response  curve. 

*  Each  data  point  includes  the  dose,  number  of  trials  at  that  dose, 

*  number  of  responses  at  that  dose,  and  the  fraction  of  trials  in 

*  which  a  response  is  observed  (=  #  responses  /  #  trials) 

*/ 

public  class  DoseResponseDataltem  implements  Comparable<DoseResponseDataItem>  { 
public  double  dose; 
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public  int  trials; 
public  int  responses; 
public  double  responseFraction; 

/* 

*  Constructors 
*/ 

DoseResponseDataltem  0  •{ 
dose  =  0.0; 
trials  =  0; 
responses  =  0; 
responseFraction  =  0; 

> 

DoseResponseDataltem  (double  dose,  int  trials,  int  responses)  { 
this. dose  =  dose; 
this. trials  =  trials; 
this. responses  =  responses; 

this .responseFraction  =  (double)  trials  /  (double)  responses 

> 

DoseResponseDataltem  (TrialDataltem  tdi)  { 
dose  =  tdi. dose; 
trials  =  1; 

responses  =  (tdi. response)  ?  1  :  0; 
responseFraction  =  (double)  responses; 

> 

/* 

*  String  representation 
*/ 

public  String  toString  ()  { 

return  String. format  ("'/.-10.2f’/,-7d7.-4d"  , 
dose,  trials,  responses); 

> 

/* 

*  Implementation  of  Comparable  interface 
*/ 

public  int  compareTo  (DoseResponseDataltem  drdi)  { 

//  order  by  ’dose’  field 

return  Double . compare  (dose,  drdi. dose); 

> 

} 


A. 4  DoseResponseData.java 

package  ng.probitf it .data; 
import  java. util.*; 

/* 
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*  Collection  of  DoseResponseDataltem  objects 
*/ 

public  class  DoseResponseData  extends  ArrayList<DoseResponseDataItem>  { 
/* 

*  Sort  the  items  in  the  collection 
*/ 

public  void  sort  ()  { 

Collections . sort  (this); 

> 

/* 

*  Factory  method  to  create  a  DoseResponseData  from 

*  a  TrialData  object 
*/ 

public  static  DoseResponseData  create  (TrialData  td)  { 

//  init  return  value 

DoseResponseData  drd  =  new  DoseResponseData  () ; 

//  use  a  HashMap,  keyed  by  the  dose  value,  to  organize  the  data 
HashMap<Double ,  DoseResponseDataItem>  map  = 

new  HashMap<Double ,  DoseResponseDataItem>  (); 

Iterator<TrialDataItem>  iter  =  td. iterator  () ; 
while  (iter .hasNext  ())  { 

TrialDataltem  tdi  =  iter. next  0 ; 

Double  key  =  new  Double  (tdi. dose); 

if  (map.containsKey  (key))  { 

//  update  #  trials,  #  responses,  and  response  fraction 
//  at  this  dose 

DoseResponseDataltem  drdi  =  map. get  (key); 
drdi .trials++; 

if  (tdi .response)  drdi .responses++; 
drdi .responseFract ion  = 

(double)  drdi . responses  /  (double)  drdi. trials; 

>  else  { 

//  entry  not  found  -  add  new  item  for  this  dose 
map. put  (key,  new  DoseResponseDataltem  (tdi)); 

> 

} 

//  insert  contents  of  map  into  DoseResponseData  and  sort 
if  (map. size  ()  >  0)  { 

Iterator<DoseResponseDataItem>  maplter  = 
map. values  (). iterator  () ; 
while  (maplter .hasNext  0)  { 
drd. add  (maplter .next  ()); 

> 

drd. sort  () ; 

> 

return  drd; 
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> 


> 


A. 5  ProbitFit.java 

package  ng.probitfit.fit; 

/* 

*  Results  of  a  Probit  Fit 

*  Includes  slope  and  intercept  of  the  fit  line, 

*  variance  of  slope  and  intercept , 

*  ybar  (weighted  avg.  of  working  probit  values)  and  its  variance, 

*  xbar  (avg.  of  loglO(dose)) , 

*  chi-squared  of  fit, 

*  #  degrees  of  freedom 

* 

*  All  these  values  are  neede  to  calculate  ED(prob)  values,  and  the 

*  associated  fiducial  limits. 

*/ 

public  class  ProbitFit  { 

public  double  a;  //  intercept 

public  double  b;  //  slope 

public  double  varA;  //  variance  of  intercept 

public  double  varB;  //  variance  of  slope 

public  double  varYbar;  //  variance  in  weighted  mean  of  probits 

public  double  chi2;  //  chi-squared  of  fit 

public  int  degFreedom;  //  degrees  of  freedom  (=  #  pts  -  2) 

public  double  xbar;  //  weighted  avg  of  loglO(dose)  (for  fid.  limits) 

public  double  ybax;  //  weighted  avg  of  probit  values  (diagnostic) 

public  int  numlterations;  //  #  iterations  required  for  fit 

public  double  Sxy;  //  diagnostic,  used  to  calc,  fit  parameters 

public  double  Syy;  //  (varB  =  1/Sxx,  varYBar  =  1/SO) 

public  boolean  infiniteSlope;  //  set  to  ’true'  if  the  slope  b  meets 

//  some  program-defined  criterion 
//  to  be  considered  ’inifinite’ 

//  This  affects  fid.  limit  calcs. 

/* 

*  Default  constructor 
*/ 

public  ProbitFit  ()  { 
a  =  0.0; 
b  =  0.0; 
varA  =0.0; 
varB  =  0.0; 
varYbar  =0.0; 
chi2  =  0.0; 
degFreedom  =  0; 
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xbar  =0.0; 
ybax  =  0.0; 
numlterations  =  0; 

Sxy  =  0.0; 

Syy  =0.0; 

inf initeSlope  =  false; 

> 

/* 

*  Copy  constructor 
*/ 

public  ProbitFit  (ProbitFit  pf)  { 
a  =  pf .a; 
b  =  pf.b; 
var A  =  pf . var  A ; 
varB  =  pf.varB; 
varYbar  =  pf . varYbar ; 
cbi2  =  pf.chi2; 
degFreedom  =  pf . degFreedom ; 
xbar  =  pf . xbar ; 
ybar  =  pf . ybar ; 

numlterations  =  pf .numlterations ; 

Sxy  =  pf.Sxy; 

Syy  =  pf - Syy ; 

inf initeSlope  =  pf. inf initeSlope; 

> 

/* 

*  String  representation,  for  debugging  purposes 
*/ 

public  String  toString  ()  { 

StringBuffer  sb  =  new  StringBuffer  ("a="); 

sb. append  (a);  sb. append  ("  "); 

sb. append  ("b=");  sb. append  (b) ;  sb. append  ("\n"); 

sb. append  ("varA=");  sb. append  (varA) ;  sb. append  ("  ") ; 

sb. append  ("varB=");  sb. append  (varB) ;  sb. append  ("\n"); 

sb. append  ("cbi2=");  sb. append  (cbi2) ;  sb. append  ("  "); 

sb. append  ("degFreedom=") ;  sb. append  (degFreedom); 

return  sb. toString  () ; 

> 

> 


A. 6  ProbitFitProcedure.java 


package  ng.probitfit.fit; 
import  java. util.*; 
import  ng. probitfit.dat a.*; 
/* 
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*  Class  that  performs  the  probit  fitting  procedured 
*/ 

public  class  ProbitFitProcedure  { 


//  default  values  for  fit  procedure  parameters 
private  static  double  DEF_TOLERANCE  =  1.0e-9; 
private  static  int  DEF_MAXITER  =  55; 

//  Slopes  greater  than  this  value  should  be  flagged  as  ’infinite’ 
private  static  double  INFINITESLOPE  =  200.0; 

//  class  fields 
private  WorkingData  wd; 

private  double  tol  =  DEF_T0LERANCE ; 

private  int  maxlter  =  DEF_MAXITER; 

private  ProbitFit  probitFit  =  new  ProbitFit  () ; 

/* 

*  Constructor  -  requires  a  DoseResponseDataSet  as  input 
*/ 

public  ProbitFitProcedure  (DoseResponseData  drd)  { 
wd  =  new  WorkingData  (drd) ; 

probitFit .degFreedom  =  wd.size  ()  -  2; 

> 

/* 

*  Accessors  to  get/set  fit  procedure  parameters 
*/ 

public  void  setTolerance  (double  tolerance)  {  tol  =  tolerance;  > 
public  double  getTolerance  ()  {  return  tol;  } 

public  void  setMaxIterations  (int  n)  {  maxlter  =  n;  > 
public  int  getMaxIterations  ()  ■(  return  maxlter;  > 

/* 

*  Return  results  of  fit 
*/ 

public  ProbitFit  getFit  ()  {  return  probitFit;  > 

/* 

*  Perform  the  fit  procedure 

*  Returns  #  iterations  if  fit  converges 

*  -1  if  max.  iterations  reached  without  convergence 
*/ 

public  int  doFit  ()  { 

//  initial  guess  for  fit  equations  is  intercept  a  =  5.0, 

//  slope  b  =  0.0,  i.e.,  all  Yp  =  5 
probitFit. a  =5.0; 
probitFit. b  =  0.0; 

//  Convergence  of  fit  process  is  determined  by  considering  the 
//  change  in  calculated  log!0(ED50)  values  between  iteration  steps 
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double  x50  =  0.0; 
double  oldX50; 


probitFit .numlterations  =  0; 
do  { 

oldX50  =  x50; 

//  update  weights,  working  probits,  etc  for  current  fit 
wd.updateForFit  (probitFit . a,  probitFit .b) ; 
iterationStep  () ; 

x50  =  (5.0  -  probitFit. a)  /  probitFit. b; 

if  (probitFit. b  >  INFINITESLOPE)  probitFit . inf initeSlope  =  true; 

++probitFit .numlterations ; 

if  (probitFit .numlterations  >=  maxlter)  { 

//  flag  no  convergence  in  prescribe  #  of  iterations  . . . 

probitFit .numlterations  =  -1; 

break;  //  ...  and  break  out  of  loop 

> 

>  while  ((Math.abs  (x50  -  oldX50)  >  tol)  &&  ! probitFit . inf initeSlope) ; 
return  probitFit .numlterations ; 

> 


/* 

*  Perform  one  step  of  the  iterative  fitting  procedure 
*/ 

private  void  iterationStep  ()  {. 

double  Sx  =  0.0;  //  sum  of  w*x 

double  Sy  =  0.0;  //  sum  of  w*ywp 

double  S  =0.0;  //  sum  of  w 

//  iterate  through  data  to  find  xbar  (avg  of  loglO(dose))  and 
//  ybar  (avg  of  working  probit  values) 
Iterator<WorkingDataItem>  iter  =  wd. iterator  () ; 
while  (iter .hasNext  ())  { 

WorkingDataltem  wdi  =  iter. next  () ; 

S  +=  wdi.w; 

Sx  +=  wdi.w  *  wdi.x; 

Sy  +=  wdi.w  *  wdi.ywp; 

> 


double  xbar  =  Sx  /  S; 
double  ybar  =  Sy  /  S; 


//  now 
//  Sxy 
double 
double 
double 
double 


iterate  to  calculate  Sxx  =  sum  (w  *  (x-xbar)"2), 

=  sum  (w*(x-xbar)*(y-ybar)) ,  Syy  =  sum  (w*(y-ybar)~2) 
Sxx  =  0.0; 

Sxy  =0.0 
Syy  =0.0 

Sb  =0.0;  //  sum  (w*(x-xbar)*ywp) 


iter  =  wd. iterator  () ;  //  reset  iterator  to  beginning  of  data 

while  (iter .hasNext  ())  { 

WorkingDataltem  wdi  =  iter. next  (); 
double  tx  =  (wdi.x  -  xbar); 
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double  ty  =  (wdi . ywp  -  ybar) ; 

Sxx  +=  wdi.w  *  tx  *  tx; 

Sxy  +=  wdi . w  *  tx  *  ty ; 

Syy  +=  wdi.w  *  ty  *  ty; 

Sb  +=  wdi . w  *  tx  *  wdi . ywp ; 

> 

//  now  calculate  intercept,  slope,  variances,  and  chi~2 

probitFit.b  =  Sb  /  Sxx; 

probitFit.a  =  ybax  -  probitFit.b  *  xbar; 

probitFit . varB  =  1.0  /  Sxx; 

probitFit . varYbar  =  1.0  /  S; 

probitFit .var A  =  1.0  /  S  +  xbar  *  xbar  /  Sxx; 

probitFit . chi2  =  Syy  -  Sxy  *  Sxy  /  Sxx; 

if  (probitFit . chi2  <  0.0) 

probitFit . chi2  =  0.0;  //  pathological  data  sets 

probitFit .xbar  =  xbar; 

probitFit . ybar  =  ybar ; 

probitFit . Sxy  =  Sxy; 

probitFit. Syy  =  Syy; 

> 

> 


A. 7  WorkingDataltem.java 

package  ng.probitfit .f it; 

import  ng .math. * ; 

import  ng.probitfit. data.*; 

/* 

*  Hold  information  from  a  single  trial  for  the  probit  fit  procedure. 

*  Holds  loglO(dose),  #  trials,  and  observed  probability  of  response. 

*  Also  holds  the  probit  value  from  the  current  fit  iteration 

*  (Y  =  a  +  b*x) ,  P(Y),  z,  weight,  and  working  probit  value 
*/ 

class  WorkingDataltem  implements  Comparable<WorkingDataItem>  { 

public  double  x;  //  loglO(dose) 

public  int  n;  //  #  trials  at  dose 

public  double  Pobs;  //  observed  #  responses  /  #  trials  at  dose 

public  double  Yp;  //  probit  value  for  current  fit  iteration 

public  double  z;  //  prob.  density  for  Yp  (from  normal  density  fn) 

public  double  Py;  //  prob.  for  Yp  (from  normal  distribution) 

public  double  w;  //  weight  for  this  point 

public  double  ywp;  //  working  probit  value 

private  static  final  double  MINP  =  10.0  *  Double ,MIN_VALUE; 

/*  constructors  */ 

public  WorkingDataltem  (double  dose,  int  n,  double  prob)  { 
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x  =  Functions . loglO  (dose); 
this.n  =  n; 

Pobs  =  prob; 

updateForFit  (5.0,  0.0); 

> 

public  WorkingDataltem  (DoseResponseDataltem  drdi)  { 
x  =  Functions . loglO  (drdi. dose); 
n  =  drdi. trials; 

Pobs  =  drdi .responseFraction; 
updateForFit  (5.0,  0.0); 


*  Calculate  probit  values,  weight,  etc.  for  a  specified 

*  intercept  a  and  slope  b 
*/ 

public  void  updateForFit  (double  a,  double  b)  { 

//  Probit  value  based  on  current  intercept  and  slope 
Yp  =  a  +  b  *  x; 

double  ypm5  =  Yp  -  5.0; 
z  =  Distributions .dnorm  (ypm5) ; 

Py  =  Distributions .pnorm  (ypm5) ; 

if  (Math.abs  (Pobs  -  Py)  <=  MINP)  { 
w  =  0.0; 
ywp  =  Yp; 

>  else  { 

w  =  (double)  n  *  z  *  z  /  (Py  *  (1.0-Py)); 
ywp  =  Yp  +  (Pobs  -  Py)  /  z; 

> 

} 

/* 

*  Comparable  interface  method 
*/ 

public  int  compareTo  (WorkingDataltem  wdi)  { 

//  order  by  x  {loglO(dose)} 
return  Double . compare  (x ,  wdi . x) ; 

> 

> 


A. 8  WorkingData.java 

package  ng . probitf it . f it ; 
import  java. util.*; 
import  ng. probitf it. data.*; 
/* 
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*  Collection  of  WorkingData,  used  to  perform  the  Probit  Fit 
*/ 

class  WorkingData  extends  ArrayList<WorkingDataItem>  { 


/* 

*  Constructor  -  creates  a  WorkingData  object  from  a 

*  DoseResponseData 
*/ 

public  WorkingData  (DoseResponseData  drd)  { 

Iterator<DoseResponseDataItem>  iter  =  drd. iterator  () ; 
while  (iter  .hasNext  0)  { 

add  (new  WorkingDataltem  (iter. next  ())); 

> 

> 

/* 

*  Sort  the  collection 
*/ 

public  void  sort  ()  { 

Collections. sort  (this); 

> 

/* 

*  Update  the  WorkingDataltem  items  in  this  collection  for  the 

*  current  iterations ’s  fit  parameters 
*/ 

public  void  updateForFit  (double  a,  double  b)  { 
Iterator<WorkingDataItem>  iter  =  iterator  0; 
while  (iter .hasNext  ())  { 

iter. next  () .updateForFit  (a,  b) ; 

> 

> 


A. 9  FitProcedureOptions.java 


package  ng.probitfit.fit; 

/* 

*  Class  to  hold  options  controlling  the  fit  procedure 
*/ 

public  class  FitProcedureOptions  { 

//  default  settings 

private  static  final  int  DEF_TOLEXPONENT  =  9; 
private  static  final  int  DEF_MAXITERATIONS  =  55; 

//  data  fields 

public  int  tolExponent;  //  fit  procedure  convergence  test  uses 

//  tol  =  10“ (-tolExponent) 

public  int  maxlterations ; 
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/* 

*  Constructors 
*/ 

public  FitProcedureOptions  ()  { 

tolExponent  =  DEF_TOLEXPONENT; 
maxlterations  =  DEF_MAXITERATIONS; 

} 

public  FitProcedureOptions  (FitProcedureOptions  fpo)  { 
tolExponent  =  fpo .tolExponent ; 
maxlterations  =  fpo. maxlterations; 

> 

> 


A.  10  EDEvaluator.java 

package  ng.probitf it .fit; 
import  ng.math. *; 

/* 

*  Class  to  calculate  ED(prob)  values  using  the  results  of  a  Probit  Fit. 

*/ 

public  class  EDEvaluator  { 

//  data  fields 

private  ProbitFit  pf ;  //  result  of  fit 

private  double  conf Level;  //  fiducial  limit  confidence  level 

private  double  heteroTestLevel;  //  test  for  heterogeneity 

private  double  h;  //  heterogeneity  factor 

private  double  t;  //  value  of  normal  or  student  distribution 

//  corresponding  to  fiducial  confidence  level 

private  double  g;  //  factor  for  fiducial  limit  calculations 

private  double  onemg;  //  (1.0  -  g) 

private  double  gdlmg;  //  g/(g-l),  factor  in  fid.  limit  calculation 
private  double  facSqrt;  //  (t*sqrt(h))/(b*(l-g)) ,  in  fid.  limit  calc, 
private  boolean  initialized;  //  flag  for  lazy  eval  of  params 
/* 

*  Constructor 
*/ 

public  EDEvaluator  (ProbitFit  probFit,  FiducialLimitOptions  options)  { 
pf  =  probFit; 

confLevel  =  options . conf idenceLevel ; 
heteroTestLevel  =  options.heterogeneityTest; 


25 


//  indicate  that  h,  t,  g,  gml  need  to  be  calculated 
initialized  =  false; 

> 

/* 

*  Read-only  accessor  for  intermediate  calculation  factors 
*/ 

public  double  getH  ()  { 
if  ( ! initialized) 
init  () ; 
return  h; 

> 

public  double  getG  ()  { 
if  ( ! initialized) 
init  ()  ; 
return  g; 


public  double  getT  ()  { 
if  (Unitialized) 
init  ()  ; 
return  t; 


/* 

*  Return  ED  and  fiducial  limits  for  requested  probablity  level 
*/ 

public  ED  getED  (double  probability)  { 
if  ([initialized) 
init  () ; 

//  get  probit  corresponding  to  probability,  recalling  that 
//  the  probit  value  is  defined  as  the  value  y  such  that 
//  pnorm  (y-5)  =  probability 

double  y  =  Distributions . qnorm  (probability)  +5.0; 

//  loglO(ED) 

double  x  =  (y  -pf .a)  /  pf . b ; 

//  logl0(f iducial  limits) 
double  xUpper  =  x; 
double  xLower  =  x; 

if  ( ! pf . inf initeSlope)  { 

//  terms  for  fiducial  limits 
double  xmxbar  =  x  -  pf . xbax ; 

double  fA; 

double  fB  =  onemg  *  pf.varYbar  +  xmxbar  *  xmxbar  *  pf.varB; 
if  (fB  <=  0.0)  { 

//  safety  check  -  not  sure  if  this  is  possible 
//  probably  need  some  really  screwed  up  data 
fB  =  0.0; 
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fA  =  0.0; 

>  else  { 

fB  =  facSqrt  *  Math.sqrt  (fB) ; 

fA  =  gdlmg  *  xmxbar;  //  g  *  (x-xbar)  /  (g-1) 

> 

//  loglO (fiducial  limits) 
xUpper  +=  (f A  +  f B) ; 
xLower  +=  (f A  -  f B) ; 

> 

ED  ed  =  new  ED  () ; 
ed.ed  =  Functions. tenTo  (x) ; 

ed.upperFL  =  Functions .tenTo  (xUpper); 
ed.lowerFL  =  Functions .tenTo  (xLower); 

return  ed; 

> 

/* 

*  Calculate  t  and  heterogeneity  factor  h 
*/ 

protected  void  init  ()  { 

double  nu  =  (double)  pf .degFreedom; 

double  px2  =  Distributions ,pchi2  (pf.ch.i2,  nu) ;  //  chi'2  prob. 

double  probLevel  =  0.5  +  conf Level/2.0; 
if  (px2  <  heteroTestLevel)  { 

//  use  heterogeneity  factor,  student  distribution 
h  =  pf.chi2  /  nu; 

t  =  Distributions. qstudent  (conf Level,  nu) ; 

>  else  { 

//  use  normal  distribution 
h  =  1.0; 

t  =  Distributions . qnorm  (probLevel); 

> 

double  temp  =  t  /  pf.b;  //  h  *  t~2  /  b~2 

g  =  h  *  temp  *  temp  *  pf.varB; 

onemg  =  1.0  -  g ; 
gdlmg  =  g  /  onemg; 

facSqrt  =  Math.sqrt  (h)  *  temp  /  onemg; 
initialized  =  true; 

> 

> 
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A. 11  ED.java 


package  ng.probitfit.fit; 

/* 

*  ED(prob)  level  and  associated  fiducial  limits 
*/ 

public  class  ED  { 

public  double  ed; 

public  double  upperFL;  //  upper  fiducial  limit 
public  double  lowerFL;  //  lower  fiducial  limit 

> 


A. 12  FiducialLimitOptions.java 

package  ng.probitfit.fit; 

/* 

*  Parameters  used  when  calculating  fiducial  limits 
*/ 

public  class  FiducialLimitOptions  { 

//  default  values  of  tbe  parameters 

private  static  final  double  DEF_CONF_LEVEL  =  0.95; 

private  static  final  double  DEF_HETERO_TEST  =  0.10; 

//  data  fields 

public  double  conf idenceLevel; 
public  double  heterogeneityTest ; 

/* 

*  Default  constructor 
*/ 

public  FiducialLimitOptions  ()  { 

conf idenceLevel  =  DEF_C0NF_LEVEL; 
heterogeneityTest  =  DEF_HETERO_TEST; 

> 

/* 

*  Copy  constructor 
*/ 

public  FiducialLimitOptions  (FiducialLimitOptions  flo)  { 
conf idenceLevel  =  flo. conf idenceLevel; 
heterogeneityTest  =  flo .heterogeneityTest; 

> 

> 


28 


A.  13  Distributions.java 


package  ng.math; 

/* 

*  Class  to  return  probability  distributions,  densities,  and  quantiles 
*/ 

public  class  Distributions  { 

//  Useful  constants 

private  static  final  double  SQRT0F2  =  Math. sqrt (2.0) ; 

private  static  final  double  SQRT0F2PI  =  Math.sqrt(2.0  *  Math. PI); 

private  static  final  double  SQRTPI  =  Math. sqrt (Math. PI) ; 

//  Class  data 

private  Functions  func  =  new  FunctionsO; 
private  static  double  tol  =  1.0e-8; 

/* 

*  Set  or  get  tolerance  for  numerical  routines 
*/ 

public  static  void  setTolerance (double  tolerance)  {  tol  =  tolerance;  } 
public  static  double  getToleranceO  {  return  tol;  } 

/* 

*  Normal  distributions 

*  pnorm(x)  =  normal  prob.  distribution 

*  dnorm(x)  =  normal  prob.  density  function 

*  qnorm(p)  =  quantile  for  normal  prob.  dist. 

*/ 

public  static  double  pnorm(double  x)  { 

return  (1.0  +  Functions .erf (x/SQRT0F2))  /  2.0; 

> 

public  static  double  dnorm(double  x)  -( 

return  Math.exp(-x*x/2.0)  /  SQRT0F2PI; 

> 

public  static  double  qnorm(double  p)  { 

//  Use  Newton-Raphson  method  to  find  the  value  x  for  which 
//  pnorm(x)  =  p 

//  Testing  indicates  #  iterations  <  10  for  values  of  interest 

//  in  the  probit  program 

double  x  =  0.0;  //  initial  guess 

double  dx; 

do  { 

dx  =  (pnorm(x)  -  p)  /  dnorm(x); 
x  -=  dx; 

>  while  (Math.abs(dx)  >  tol); 
return  x; 

> 
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/* 

*  Chi-square  distribution 

*  Probability  of  observing  a  larger  value  of  chi-squared  than  the 

*  input  parameter  chi2  for  the  nu  degrees  of  freedom 

*  pchi2  large  ->  prob  of  observing  a  larger  chi-squared  is  high 

*  ->  good  fit 
*/ 

public  static  double  pchi2 (double  chi2,  double  nu)  { 
return  1.0  -  Functions . gammap (nu/2 . 0 ,  chi2/2.0); 

> 

/* 

*  Student’s  distribution 

*  pstudent  -  prob.  function 

*  qstudent  -  quantile 
*/ 

public  static  double  pstudent (double  t,  double  nu)  { 
double  x  =  nu  /  (nu  +  t*t) ; 

return  1.0  -  Functions .betaine (nu/2. 0,  0.5,  x) ; 

> 

public  static  double  qstudent (double  p,  double  nu)  { 

//  pre-calculate  some  common  factors  for  the  derivative 

double  nupld2  =  (nu  +  1.0)/2.0; 

double  nud2  =  nu/2.0; 

double  facl  =  Math.log(2.0/SQRTPI) 

+  Functions. gammaln(nupld2) 

-  Functions. gammaln(nud2) 

+  nud2*Math.log(nu) ; 

double  q  =  0.0;  //  initial  guess 
double  dq; 
do  { 

double  fac2  =  facl  -  nupld2*Math. log(nu  +  q*q) ; 
dq  =  (pstudent (q,nu)  -  p) /Math. exp (f ac2) ; 
q  -=  dq; 

}  while  (Math . abs (dq)  >  tol) ; 
return  q; 

> 


A. 14  Functions.java 

package  ng.math; 

/* 

*  Class  provides  special  functions 

* 

*  Includes: 

*  erf(x)  Error  function 
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*  loglO(x)  Log  based  10 

*  tenTo(x)  10'x 

*  gammaln(x)  log  of  Gamma  function 

*  gammap(a,x)  incomplete  gamma  function 

*  betaine (a, b,x)  incomplete  beta  function 
*/ 

public  class  Functions  { 

//  Useful  constants 

private  static  final  double  SQRTPI  =  Matb. sqrt (Math. PI) ; 
private  static  final  double  SQRT2PI  =  Math.sqrt(2.0*Math.PI) ; 
private  static  final  double  L0G0F10  =  Math. log ( 10. 0) ; 
private  static  final  double  TINY  =  1.0e-30; 

//  Class  data 

private  static  double  tol  =  1.0e-14;  //  tolerance  for  numerical  routines 

/* 

*  Set  or  get  tolerance  for  numerical  calculations 
*/ 

public  static  void  setTolerance (double  tolerance)  {  tol  =  tolerance;  } 
public  static  double  getTolerance 0  {  return  tol;  } 

/* 

*  erf(x)  -  Error  function 
*/ 

public  static  double  erf (double  x)  { 
double  sign  =  1.0; 
if  (x  <  0.0)  { 

x  =  Math .  abs  (x)  ; 
sign  =  -1.0; 

> 

//  Testing  indicates  series  solution  requires  fewer  iteration  for 
//  x  <  2.5,  while  continued  fraction  converges  faster  for  x  >  2.5 
if  (x  <  2.5)  { 

return  sign  *  erf Series (x) ; 

>  else  { 

return  sign  *  (1.0  -  erfcContFrac(x)) ; 

> 

> 

//  Series  solution  for  error  func  -  called  by  erf(x) 
private  static  double  erf Series (double  x)  { 

//  Series  expansion  from  Abramowitz  and  Stegun 
int  k  =  0; 

double  an  =  1.0;  //  series  term  coeff. 

double  sum  =  1.0;  //  value  from  1st  term  of  series 

double  oldsum; 

//  Not  checking  #  iterations  here  -  the  erf(x)  routine  only  calls  this 
//  method  for  values  of  x  <  2.5,  a  region  for  which  at  max  30-40 
//  iterations  are  required 
do  { 
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oldsum  =  sum; 

++k; 

an  *=  -x*x; 

an  /=  (double)  k; 

sum  +=  an  /  (2.0  *  k  +1.0); 

>  while  (Math. abs (sum  -  oldsum)  >  tol) ; 

return  2.0  *  x  *  sum  /  SQRTPI; 

> 

//  Continued  fraction  solution  for  complimentary  error  function  - 
//  called  by  erf (x) 

private  static  double  erf cContFrac (double  x)  { 

//  CF  expression  from  Abramowitz  and  Stegun 
//  Lentz’s  method  to  evalute  CF  base  on  material  from 
//  Numerical  recipes  in  C 

double  c  =  x  +  1.0/TINY;  //start  from  bO  =  0 
double  d  =  1.0/x; 
double  delta  =  c  *  d; 

double  f  =  TINY  *  delta;  //  This  takes  care  of  1st  iteration 
int  i  =  1; 

double  a  =  1.0;  //  all  bn  =  x  for  n  >  0 
do  { 

++i; 

a  =  (double)  (i-1)  /  2.0; 

//  recursion  relations  for  cn  and  dn 
d  =  x  +  a*d; 

if  (d  ==  0.0)  d  =  TINY;  //  Dangerous  check 
d  =  1.0/d; 
c  =  x  +  a/c; 

if  (c  ==  0.0)  c  =  TINY;  //  Dangerous  check 
delta  =  c  *  d; 
f  =  f  *  delta; 

}  while  (Math. abs (delta- 1 .0)  >  tol); 
return  f  *  Math. exp (-x*x) /SQRTPI; 

> 

/* 

*  loglO(x)  -  log  based  10  of  x 

*  tenTo(x)  -  10~x 
*/ 

public  static  double  loglO (double  x)  { 
return  Math. log (x)  /  L0G0F10; 

> 

public  static  double  tenTo (double  x)  { 
return  Math. exp (x  *  L0G0F10) ; 

> 

/* 

*  Natural  log  of  gamma  functions 
*/ 
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private  static  final  double  gammaln_coeff []  =  { 

76.18009172947146, 

-86.50532032941677, 

24.01409824083091, 

-1.231739572450155, 

0. 1208650973866179e-2, 

-0 . 5395239384953e-5 

>; 

public  static  double  gammaln (double  x)  { 

//  Routine  adapted  from  Numerical  Recipes  in  C 

double  temp  =  x  +  5.5; 

temp  -=  (x+0.5)*Math.log(temp) ; 

double  y  =  x; 

double  sum  =  1.000000000190015; 
for  (int  i  =  0;  i  <  6;  i++)  { 

++y; 

sum  +=  gammaln_coeff [i] /y ; 

> 

return  -temp  +  Math. log (SQRT2PI  *  sum  /  x) ; 

> 

/* 

*  Incomplete  gamma  functions  P(a,x) 

*/ 

public  static  double  gammap (double  a,  double  x)  { 

//  Routine  adapted  from  Numerical  Recipes  in  C 
if  (x  <  0.0)  { 

throw  new  IllegalArgumentException( "gammap:  Illegal  arg,  x="  +  x) ; 

} 

if  (a  <  0.0)  { 

throw  new  IllegalArgumentExceptionO'gammap:  Illegal  arg,  a="  +  a); 

} 

if  (x  <  a  +  1.0)  { 

return  gamSeries(a.x) ; 

>  else  { 

return  1.0  -  gamContFrac (a,x) ; 

> 

> 

//  Support  functions  for  gammap 

private  static  double  gamSeriesTolerance  =  1.0e-10; 
private  static  double  gamSeries (double  a,  double  x)  { 

//  Series  solution  for  P(a,x) 

double  t  =  1.0/a; 

double  sum  =  t; 

double  oldsum; 

double  aa  =  a; 

do  ■( 

oldsum  =  sum; 

++aa; 

t  *=  (x/aa) ; 
sum  +=  t; 

}  while  (Math.abs(sum  -  oldsum)  >  gamSeriesTolerance); 
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> 


return  sum  *  Math.exp(-x  +  a*Math.log(x)  -  gammaln(a)); 


private  static  int  GAMCF_MAXITER  =  100; 

private  static  double  gamContFrac (double  a,  double  x)  { 

//  Constant  fraction  for  1.0  -  P(a,x) 

//  Routine  adapted  from  Numerical  Recipes  in  C 

double  b=x+1.0-a; 

double  c  =  1.0/TINY; 

double  d  =  1.0/b; 

double  h  =  d; 

double  delta; 

int  niter  =  0; 
do  { 

++nlter ; 

double  an  =  -niter  *  (niter  -  a) ; 
b  +=  2.0; 
d  =  an*d  +  b; 

if  (Math . abs (d)  <  TINY)  d  =  TINY; 
c  =  b  +  an/c; 

if  (Math. abs (c)  <  TINY)  c  =  TINY; 
d  =  1.0/d; 
delta  =  d*c; 
h  *=  delta; 

if  (niter  >=  GAMCF_MAXITER)  { 

throw  new  RuntimeExceptionC'gamConstFrac :  No  convergence"); 

> 

>  while  (Math. abs (delta  -  1.0)  >  tol) ; 

return  Math.exp(-x  +  a*Math. log(x)  -  gammaln(a))*h; 

> 

/* 

*  Incomplete  beta  function  Ix(a,b) 

*/ 

public  static  double  betaine (double  a,  double  b,  double  x)  { 

//  Routine  from  Numerical  Recipes  in  C 
if  (x  <  0.0  I  I  x  >  1.0)  { 

throw  new  IllegalArgumentExceptionO'betainc:  Illegal  arg  x="  +  x) ; 

} 

double  bt; 

if  (x  ==  0.0  II  x  ==  1.0)  { 
bt  =  0.0; 

}  else  { 

bt  =  gammaln ( a+b) -gammaln (a) -gammaln (b) 

+  a*Math. log(x)+b*Math. log(l . 0-x) ; 
bt  =  Math. exp (bt ) ; 

> 

if  (x  <  (a+1.0)/(a+b+2.0))  { 

//  Use  continued  fraction  directly 
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return  bt  *  betaContFrac(a,b,x)/a; 

}  else  { 

//  Use  continued  fraction  after  making  symmetry  transformation 
return  1.0  -  bt  *  betaContFrac (b , a , 1 . 0-x) /b ; 

> 

> 

//  Continued  fraction  for  incomplete  beta 
private  static  int  BETACF_MAXITER  =  100; 

private  static  double  betaContFrac (double  a,  double  b,  double  x)  { 

//  Routine  from  Numerical  Recipes  in  C 
double  qab  =  a  +  b; 
double  qap  =  a  +  1.0; 
double  qam  =  a  -  1.0; 

double  c  =  1.0; 

double  d  =  1 . 0  -  qab*x/qap ; 

if  (Math . abs (d)  <  TINY)  d  =  TINY; 

d  =  1.0/d; 

double  h  =  d; 

int  niter  =  0; 
double  delta; 
do  { 

++nlter; 

double  m2  =  2.0  *  niter; 

double  aa  =  niter* (b-nlter) *x/ ( (qam+m2) * (a+m2) ) ; 

//  One  step  (even  one)  of  the  recurrence 
d  =  1 . 0  +  aa*d ; 

if  (Math. abs (d)  <  TINY)  d  =  TINY; 
d  =  1.0/d; 
c  =  1.0  +  aa/c; 

if  (Math. abs (c)  <  TINY)  c  =  TINY; 
h  *=  d*c; 

//  One  step  (odd  one)  of  the  recurrence 

aa  =  -(a+nlter)*(qab+nlter)*x/((a+m2)*(qap+m2)) ; 

d  =  1.0  +  aa*d; 

if  (Math. abs (d)  <  TINY)  d  =  TINY; 
d  =  1.0/d; 
c  =  1 . 0  +  aa/  c ; 

if  (Math. abs (c)  <  TINY)  c  =  TINY; 
delta  =  d*c; 
h  *=  delta; 

if  (niter  >=  BETACF_MAXITER)  { 

throw  new  RuntimeExcept ion ("betaContFrac:  No  convergence") 

> 

}  while  (Math. abs (delta  -  1.0)  >  tol) ; 
return  h; 


> 


35 


